% FFT_Operator obtainment: L,DL,PL
function [ts,L,DL,PL] = FFT_Operator(NP,ND,NH,v,w)

%% ------------------L,DL field start------------------------------
	T=2*pi*v/w;
	DT=T/NP;
	ts=(0:DT:T)';                                                           		% the length of ts is NP ts=(0:DT:T-DT)';
	L=zeros((NP+1)*ND,(2*NH+1)*ND);
	DL=zeros((NP+1)*ND,(2*NH+1)*ND);                                               		% dL/dt
	L(1:(NP+1)*ND,1:ND)=kron(eye(ND),ones(NP+1,1).*1/sqrt(2));
	DL(1:(NP+1)*ND,1:ND)=kron(eye(ND),zeros(NP+1,1));
	for jj=1:NH
		L(1:(NP+1)*ND,(2*jj-1)*ND+1:2*jj*ND)=kron(eye(ND),sin(jj*w/v*ts));
		L(1:(NP+1)*ND,2*jj*ND+1:(2*jj+1)*ND)=kron(eye(ND),cos(jj*w/v*ts));
		DL(1:(NP+1)*ND,(2*jj-1)*ND+1:2*jj*ND)=kron(eye(ND),(jj*w/v)*cos(jj*w/v*ts));
		DL(1:(NP+1)*ND,2*jj*ND+1:(2*jj+1)*ND)=kron(eye(ND),-(jj*w/v)*sin(jj*w/v*ts));   
	end
% ------------------L,DL field end------------------------------
	PL = pinv(L);
end
